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Abstract 

PACS number 

The method is an extension to negative energies of a spectral integral equation method to 

Q ■ solve the Schroedinger equation, developed previously for scattering applications. One important 

Q , innovation is a re-scaling procedure in order to compensate for the exponential behaviour of the 

negative energy Green's function. Another is the need to find approximate energy eigenvalues, to 

Q_i' serve as starting values for a subsequent iteration procedure. In order to illustrate the new method, 

the binding energy of the He-He dimer is calculated, using the He-He TTY potential. In view of 

, the small value of the binding energy, the wave function has to be calculated out to a distance 

' of 3000 a.u. Two hundred mesh points were sufficient to obtain an accuracy of three significant 

■ figures for the binding energy, and with 320 mesh points the accuracy increased to six significant 

^ . 

O . figures. An application to a potential with two wells, separated by a barrier, is also made. 
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I. INTRODUCTION 



The much used differential Schrodinger equation is normally solved by means of a fi- 
nite difference method, such as Numerov or Runge-Kutta, while the equivalent integral 
Lippmann-Schwinger (LS) equation is rarely solved. The reason, of course, is that the for- 
mer is easier to implement than the latter. However, a good method for solving the LS 
equation has recently been developed 1], and applications to various atomic systems have 
been presented This method, denoted as S-IEM (for Spectral Integral Equation 

Method), expands the unknown solution into Chebyshev polynomials, and obtains equa- 
tions for the respective coefficients. The expansion is called "spectral", because it converges 
very rapidly, and hence is economical in the number of meshpoints required in order to 
attain a prescribed accuracy. A basic and simple description of the method has now been 
published 0], and a MATLAB implementation is also included. However, the applications 
described so far refer to positive energies, i.e., to scattering situations, while an example for 
negative energies, i.e., bound states, has up to now not been provided. 

Since the solution of many quantum-mechanical problems requires the availability of a ba- 
sis of discrete negative energy eigenfunctions (or bound states), or of positive energy Sturm- 
Liouville eigenfunctions, the S-IEM has now been adapted to also obtain eigenfunctions 
and eigenvalues. Since there are situations where the commonly known eigenvalue-finding 
methods do not work well, we present here a short description of our method, in the hope 
that it will be useful for the physics student/teacher community. 

An illustration of the method for the case of the bound state of the He-He atomic dimer 
is presented. This is an interesting case since the binding energy is very small, 1.3 mK 
or 1.1 X lO^'^eV , and the corresponding wave function extends out to large distances, 
between 1000 to 3000 atomic units, depending on the accuracy required. Hence a method is 
desirable that maintains accuracy out to large distances, and that can find small eigenvalues. 
A commonly used method to obtain eigenvalues consists in discretizing the Schrodinger 
differential operator into a matrix form, and then numerically obtaining the eigenvalues of 
this matrix. This procedure gives good accuracy for the low-lying (most bound) eigenvalues, 
while the least bound eigenvalues become inaccurate. The method described here does not 
suffer from this difficulty since it finds each eigenvalue of the integral equation iteratively 
It also provides a good search procedure for finding initial values of the eigenvalue, required 
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to start the iteration. 



II. THE FORMALISM. 

For negative energy eigenvalues the differential equation to be solved is 

where f is the radial distance in units of length, V and E are the potential energy and the 
(negative) energy in units of energy, respectively. This is the radial equation for the partial 
wave of angular momentum 0. For atomic physics applications this equation can be written 
in the dimensionless form 

-'^ + iV + .^)^ = (2) 

where r = f/ao is the relative distance in units of Bohr, and V and are given in atomic 
energy units. The LS eigenvalue equation that is the equivalent to Eq. is 

^(r) = [ g{r,r')V{r')i;{r')dr' (3) 
Jo 

where, as is well known, the Green's function Q{r, r') for negative energies —E = {2M/h'^) 
is given by 

g{r,r') = --E{r^)G{r'^) (4) 
r< and r> being the lesser and larger values of r and r', respectively, and 



E{r) = sinh(Kr), G{r) = exp(— /tr). (5) 

The Eq. is a Fredholm integral eigenvalue equation of the first kind. Unless the 
wave number k has a correct value, the solution does not satisfy the boundary condition 
that iplr) decay exponentially at large distances. As shown by Hartree many years ago, 
a method of finding a correct value of k, is to start with an initial guess Kg for k, divide 
the corresponding (wrong) wave function into an "out" and and "in" part, and match the 
two at an intermediary point Tm- The out part tpo is obtained by integrating (j2)) from the 
origin to an intermediate radial distance Tm, and ipj is the result of integrating Q from the 
upper limit of the radial range T inward to Tm- For the present application the integration 
method is based on the S-IEM, described in Appendix 1. The function iJjq is renormalized 



so as to be equal to ipi aX r = Tm and its value at r = Tm is denoted as ijjM- The derivatives 
with respect to r at r = Tm are calculated, as described in Appendix 1, and are denoted as 
tpQ and ip'j, respectively. The new value of the wave number k^+i is given in terms of these 
quantities as 

Ks+i = Ks- {Iter)s (6) 

where 

{Iter). = J_ ^m{^o-^'i)m 

Equations © and (|7j) can be derived by first writing Q for the exact wave function ipE 
(using for and © for the approximate wave function tpA = (V'o or ipi), multiplying 
each equation by the other wave function, integrating over r, and subtracting one from 
another. When Kqo is replaced by Hg+i and i/'e is replaced by ipA then equations (jH} and (jZ} 
result. 



III. THE SPECTRAL METHOD 



The S-IEM procedure to evaluate ifjo and ipi is as follows. First the whole radial interval 

< r < T is divided into m partitions, with the i-th partition defined as < r < ti, 

1 = 1, 2, ■ ■ - m. For notational convenience we denote the i-th partition simply as i. In each 
partition i two independent functions yi{r) and Zi{r) are obtained by solving the integral 
equations 

y,{r) = j ' g{r, r')V{r')y,{r')dr' + /,(r) (8) 

and 

z^{r)= r g{r,r')Vir')z^ir')dr' + g,{r). (9) 
J ti-1 

Here fi and Qi are scaled forms of the functions F and G defined above on the interval i, 

fi{r) = sinh(Kr) x Ei, gi{r) = exp(-Kr)/Ei, (10) 

and the scaling factor Ei in each partition i is given by 

Ei = exp{-Kti). (11) 

Such scaling factors are needed in order to prevent the unsealed functions sinh(Kr) and 
exp(— Kr), and the corresponding functions Yi and Zi to become too disparate at large dis- 
tances, which in turn would result in a loss of accuracy. Apart from these scaling operations, 
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the calculation of functions i/i and Zi by means of expansions into Chebychev polynomials, 
as well as the determination of the size of the partition i in terms of the tolerance parameter 
e is very similar to the calculation of the functions 1^ and described in Ref. (P]). The 
number of Chebychev polynomials in each partition is normally taken as = 16. The 
equations (jH)) and Q are Fredholm integral equation of the 2nd kind, and hence are much 
easier to solve than the Fredholm equations of the first kind. 
The global wave function ip is given in each partition by 

tpir) = aiyi{r) + hiZi{r). (12) 

In order to obtain the coefficients and hi for each partitions i one proceeds similarly to 
"Method B" described in Ref. (0]), that relates these coefficients from one partition to 
those in a neighboring partition. That relation is 



Ei/Eij^i 
















= li 




Ei+i/Ei 




bi+1 




h 



(13) 



where the elements of the 2x2 matrices uj and 7 are given in terms of overlap integrals 
{fy)i^ {fz)i, {gy),, {gz)^, of the type {fy). = J^^^^ fi{r)V{r)yi{r)dr, as is described in 
further detail in the Appendix 1. This relation enables one to march outward by obtaining 
ao,i+i and bo,i+i in terms of ao,i and bo,i, and inward by obtaining a/.j and bj^i in terms of 
a/^i+i and 6/^j+i.The integration outward is started at the innermost partition i = 1 with 
(^0,1 = ^/Ei, and the integration inwards is started at the outermost partition (ending at 
T), for which the coefficients and bm are given as and Em, respectively. The values of 
the functions / and O and their derivatives at the inner matching point Tm, as well as the 
integrals J^'^' ippdr + ipjdr, required for evaluating Iter in Eq. ((Zj), can be obtained in 
terms of the overlap integrals described above, as is described in Appendix 1. The iteration 
for the final value of k proceeds until the value of Iter is smaller than a prescribed tolerance. 
The important question of how to find an initial value of k is described in the next section. 

IV. SEARCH FOR THE INITIAL VALUES OF k 

Since the present method does not obtain all the values of the energy as the eigenvalues 
of one big matrix, but rather obtains iteratively one selected eigenvalue at a time, it is 
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necessary to have a reliable algorithm for finding the appropriate starting values kq for the 
iteration procedure. 

The present search method is based on Eq. (jl2|) . according to which the solution ip in 
a given partition i is made up of two parts, yi{r) and Zi{r). In the radial regions where 
the potential is small compared to the energy, i.e., in the "far" region beyond the outer 
turning point, the functions yi{r) and Zi{r) are nearly equal to the driving terms / and g 
of the respective integral equations (jH)) and 0. Hence, for negative energies, according to 
Eqs. (fTUI) . in the "far" region yi{r) has an exponentially increasing behavior, while Zi{r) is 
exponentially decreasing. For the correct bound state energy eigenvalue the solution ip has 
to decrease exponentially at large distances, and hence the coefficient Oj in Eq. (|12j) has 
to be zero for the last partition i = m. Hence, as a function of k the coefficient goes 
through zero at a value of n equal to one of the the bound state energies. 

Based on the above considerations, the search procedure for the initial value kq is as 
follows: A convenient grid of equispaced Kg values is constructed, s = 1, 2, ... and for each Kg 
the integration "outward" for the wave function is carried out to Tm ^ T , but Iter is not 
calculated. The value of Tm is selected such that the potential V is less than the expected 
binding energy. The values of the coefficient ao,iM ^^^t partition are recorded, 

and the values of Kg for which ao.j^/ changes sign are the desired starting values kq for the 
iteration procedure. The numerical example, given in the sections describing the calculation 
of the He — He bound state, shows that this search method is very reliable. 



A. The Numerical Code 



The code was written in MATLAB, and is available from the authors both in MATLAB 
and in FORTRAN versions. The code that performs the iterations is denoted as Iter_neg-k, 
and the search code for finding the starting values kq is denoted as Searchah_negJi. The 
subroutines for both codes are the same. The validity of the code was tested by comparing 
the resulting binding energy with a non-iterative spectral algorithm that obtains the eigen- 
values of a matrix. The potential used for this comparison was an analytical approximation 
to the He — He potential TTY [3], described in the next section. The comparison algorithm 
expands the wave function from Rstart to T (no partitions) in terms of scaled Legendre 
polynomials up to order A^. The operator —dP/dr'^ + \^ is discretized into a matrix at zeros 
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of the Legendre polynomial of order + 1 . The boundary conditions that the wave function 
vanishes at both Rstart and at T are incorporated into the matrix, and the eigenvalues of 
the matrix are calculated. The agreement between the two codes for the binding energy was 
good to 6 significant figures. 

In the test-calculation for the He — He dimer binding energy described below, the con- 
vergence rate of the iterations, the stability with respect to the value of a repulsive core 
cut-off parameter, and also the number of mesh-points required for a given input value of 
the tolerance parameter will be examined. A bench-mark calculation of the dimer binding 
energy is also provided for students that would like to compare their method of calculation 
to ours. In these calculations the TTY potential is replaced by an analytical approximation 
that is easier to implement. 



V. APPLICATION TO THE He - He DIMER 



The He-He dimer is an interesting molecule, because, being so weakly bound, it is the 
largest two-atom molecule known. The He-He interaction, although weak, does influence 
properties such as the superfluidity of bulk He II, of He clusters, the formation of Bose- 
Einstein condensates, and the calculation of the He trimer. In 1982 Stwalley et al |J] were 
the first to conjecture the existence of a He- He dimer. The first experimental indication of the 
dimer's existence was found in 1993 Q], and since 1994 it was explored by means of a series 
of beautiful diffraction experiments. Through these diffraction experiments not only has the 
existence of the dimer, but also that of the trimer, been unequivocally demonstrated and an 
indication of the spatial extent of these molecules has also been obtained [3], 0], Various 
precise calculations of the He- He interaction have subsequently been performed [10] and the 
corresponding theoretical binding energies of the dimer (close to 1.3 niK ^ 1.1 x 10~'^eV, see 
Table 1 in Ref. |ll|'l and the trimer (the ground state of the trimer is close to 126mK, see 
for instance Ref. 12]) agree with experiment to within the experimental uncertainty. The 
wave function of the He dimer or trimer extends out to large distances (several thousand 
atomic units), the binding energy is very weak, and the transition from the region of the 
large repulsive core to the weak attractive potential valley is very abrupt. For these reasons 
the dimer (or trimer) calculations involving He atoms require good numerical accuracy, and 
therefore was chosen as a test case for our new algorithm. 
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FIG. 1: The "TTY" He-He potential given by Tang, Toennies, and Yiu [ij, and the fit FIT 4 , as 
a function of distance. 



The transition from Eq. to the dimensionless Eq. ((21), is accomplished by transforming 
the potential and the energy into dimensionless quantities as follows 



where Q is a normalization constant, defined in Appendix 2. For the case of two colliding 
He atoms interacting via the TTY potential we take the mass of the He atom as given in 
Ref. for which the value of Q is 7295.8356. For the calculations involving our analytical 
fits to the TTY potential, we take for Q the value 7296.3. 



The TTY potential |l3|, and one analytic fit, are shown in Fig. (^. The repulsive core 
goes out to about 5 and the subsequent attractive valley reaches its maximum depth of 
3.5 X lQ~^au (approximately lO^^eV^) near r ~ 5.6 Oq. This attractive potential valley then 
decays slowly over toge distances approximately like The corresponding energy of the 
bound state is ~ —10 '^ev [13]. In the units defined in (j2I) the potential valley has a depth 
of 0.26 and the binding energy has the value of 3.04 x 10^^. The bound state wave function 
peaks near r = 10 ao and decays slowly from there on. The outer turning point occurs near 



V = QV 



(14) 



= -QE 



(15) 
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B.E.{mK) 


< r > (nm) 


Present 


1.31461 


5.1607 


Ref. 


1.30962 




Ref. \13] 


1.316 




Experiment [8] 


0.9 - 1.4 


5.2 ±0.4 



TABLE I: Comparison of the He-He Binding Energies obtained by various authors. 

30 ao; the value of the wave function at r = 2500 is ~ 10~^, and at 3000 ao it is ~ 6 x 10~^. 
The quantity r x ip"^ has its maximum beyond the turning point near r = 100 Oq, and the 
average radial separation (r) = ip"^ r dr is close to ~ 97 Oq. 

A. Results for the TTY Potential 

The TTY He-He potential is calculated by means of Fortran code provided by Franco 
Gianturco and modified at hoc for small distances (less than 1 oq) so that it maintains 
the repulsive core nature. The potential is "cut off" at a distance Rcut so that for r < Rcut, 
V{r) = V{Rcut)- The S-IEM calculation starts at r = and extends to T = 3,000 ao. The 
intermediary matching point is Tm = 7 a^. The dependence of the eigenvalue on Rcut, and 
the rate of convergence of the iterations, are described in Appendix 3. Our choice for the 
value of Rcut = 2.5 ao, of T = 3, 000, and of the tolerance parameter e = 10^^^ is such that 
the numerical stability of our results is better than 12 significant figures. 

Our value for the binding energy is compared with that of other calculations in Table 
ni. The comparison shows good agreement of our result with the literature. The difference 
between our S-IEM result and that of Ref. could well be due to a slightly different 
choice of the parameters that determine TTY. 

B. Numerical Properties of the S-IEM. 

In order to examine the nature of the partition distribution and the resulting accuracy as a 
function of the tolerance parameter e and also in order to provide a bench-mark calculation, 
the TTY potential was replaced by an analytical approximation defined in the equation 
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parameters 


FIT 3 


FIT 4 


Pi 


-3.4401e-5 


-2.930 e-5 


P2 


5.606 


5.590 


P3 


0.8695 


0.8511 


PA 


7.657 


7.5892 


P5 


1.750 


0.95608 


P6 


0.6784 


0.89098 



TABLE II: Parameters for two analytic fits to the TTY potential 

below. 

V{r) = pi exp(-(r - P2)M) x [2 - exp(-(r - P2)/P3)] 

- P6 X (r-^-«°^)/ {1 + exp [-p5(r - ^4)]} ■ (16) 

The parameters pi to for two fits, denoted FIT 3 and FIT 4, are given in Table 2. The 
resulting potential is in atomic units; r,p2,P3 and p^ are in units of oq; P5 is in units of 
and pi and pq are in atomic energy units. For all calculations involving these analytical fits, 
Q defined in Eq. (fT^ . has the value 7296.3. Fits 3 (4) produce a more (less) repulsive core 
than TTY, and is more (less) attractive in the region of the potential minimum. 

Our algorithm automatically chooses the size of the partitions such that the error in 
the functions calculated in each partition does not exceed the tolerance parameter e. At 
small distances the density of partitions is very high, but beyond 500 ag the size of the 
partitions increase to about 440 oq. In the region near the repulsive core the partitions are 
approximately 0.5 Oq wide, but there is a region in the vicinity of Rcut where they crowd 
together much more. The latter is illustrated in Fig. (0) for FIT 4, for Rcut = 2.5 oq, for 
various values of the tolerance parameter e. In Table UTTl the corresponding accuracy of the 
binding energy is displayed, for the case that the He — He potential energy is given by Fit 
4. It is noteworthy that the number of reliable significant figures in k tracks faithfully the 
value of the tolerance parameter, as is shown in Table IIIIl 
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FIG. 2: Partition distribution in the radial region up to 4 oq for three different values of the 
tolerance parameter. The value of the latter is listed in the legend and in Table HTll The potential 
is given by FIT 4, described in the text and the value of Rcut = 2.5 ao. The total number of 
partitions for each case is given by the numbers near the curves. 



Tol 


K X 10^ (ao)"^ 


M 


No. of Adeshpts 


10-12 


5.0817542 


47 


652 


10-6 


5.0817461 


19 


275 


10-^ 


5.0776 


13 


208 



TABLE III: Accuracy of the wave number (it is the square root of the binding energy) as a function 
of the tolerance parameter. The total number of partitions for each case is denoted by M. The 
corresponding partition distributions are displayed in Fig. 4 

C. The Search for the Starting Values kq. 

An example of the search procedure is given in Table IIVI for a potential given by Fit 
4, Eq. (fT^ . multiplied by the factor A = 20. The mesh of k values starts at k = —2 and 
proceeds by steps of Ak = 0.05 until k = —0.05 (all in units of OqI). The mesh values of 
K for which the coefficient a^n of ym{f^) changes sign are shown in the first column of Table 
IIV| and the corresponding iterated value of n is shown in the third column. The value of 
Tm = 80 A.M., and iter = 10^^. The MATLAB computing time required for carrying out 
the 40 mesh search calculations is 3.8 s on a 2 GHz PC; the approximately 7 iterations 
required for obtaining the more precise values of each k shown in the third column take 
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K — Mesh 


sign of am 


K — Iterated 


# of nodes 


1.70 


+ ^ - 


1.7028 





0.70 




0.7273 


1 


0.05 


+ ^ - 


0.0561 


2 



TABLE IV: Search of the wave number eigenvalues for a He-He Fit 4 potential multiplied by 20 



FIG. 3: The eigenvalues of k as a function of the strength parameter A of the Fit 4 Potential, 
approximately 1 s. of computer time.. 

By repeating the same procedure for different values of A, one can trace the k eigenvalues 
down to A = 1. The result, displayed in Fig. (jH)), shows that the values of n depend nearly 
linearly on the value of A. Futher searches with values of A slightly less than unity showed 
that the code was able to find an energy that is approximately 25 times less bound than the 
result for the TTY potential. 

In order to provide a benchmark calculation, the values of n obtained with the potential 
of FIT 4 are listed in Table IVl The value of Q is 7296.3, the values of Tm and T are 7 and 
3, 000 ao, respectively, and A = 1.0 
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Rcut (q^o) 




2.5 


5.08175419E-3 


3.0 


5.08176556E-3 


3.5 


5.10608688E-3 



TABLE V: The values of the bound state wave number for potential of FIT 4, for various values 
of the cut-off radius . Additional information is given in the text 

VI. APPLICATION TO A DOUBLE WELL POTENTIAL. 

The case for which the potential has two (or more) wells separated by one (or more) 
barriers offers another test for the reliability and accuracy of a numerical procedure for 
obtaining eigenvalues of the Schrodinger equation. The reason is that the energy eigenvalues 
are split by a small amount, corresponding to the situation in which the wave function located 
in one of the wells has either the same or the opposite sign of the wave function located 
in the adjoining well. The larger the barrier, the smaller is the difference AE between the 
two energies, and the larger are the demands on the numerical procedure. An interesting 
relaxation method for finding energy eigenvalues contained in a prescribed interval has been 
described in Ref . |3| • The double well potential, in the units of Eq. (j2)) is 

V = -Ax^ + x\ -Tm<x< T^. (17) 

The value of AE for the difference between the two lowest eigenvalues were calculated here by 
using the S-IEM method described in this paper, denoted as (Ai?)/£;Mand also by a matrix 
eigenvalue method, denoted as {AE)l. This method discretizes the Schrodinger operator on 
the left hand side of Eq. ^ 

(-|^ + ^)^ = ^^ (18) 
at the zeros of a Legendre Polynomial of order til, and then finds all the eigenvalues of the 
corresponding matrix using the standard QR algorithm. The comparison of the results for 
the three largest values of A is shown in Table EU where {AE)rei denotes the result obtained 
in Ref. |3| For the S-IEM the value of the tolerance parameter was e = 10~^^, and the 
corresponding accuracy was sufficient to obtain the results shown in the Table IVTl However, 
it can be seen that the relaxation method is more accurate than the S-IEM method. The 
difference between the results in Table |VT1 could well be due to differences in the choice of 
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A 








10 


3.02E-5 


2.98185E-5 


2.9821E-5 


12 


3.53E-7 


3.508E-7 


3.5093E-7 


15 


2E-10 


2E-10 


1.9499E-10 



TABLE VI: Comparison between three different metliods of calculating energy eigenvalues. The 
table shows the difference between the two lowest eigenvalues of the double well potential defined 
in this section 

the value of Tm- For the {AE)l result, the value of Tm was varied between 6 and 9 units of 
length, and the number of Legendre polynomials rii was varied between 200 and 700. The 
numerical stability of the QR algorithm is well documented in the numerical linear algebra 
literature. The convergence of the Legendre discretization of the Schrodinger operator using 
finite series expansions in orthogonal polynomials, such as Legendre, Chebyshev and others, 
is also well understood, as discussed for example in Ref 



VII. SUMMARY AND CONCLUSIONS. 

An integral equation method (S-IEM) ^ for solving the Schrodinger equation for positive 
energies has been extended to negative bound-state energy eigenvalues. Our new algorithm 
is in principle very similar to an iterative method given by Hartree in 1930, in that it guesses 
a binding energy, integrates the Schrodinger Equation inwards to an intermediary matching 
point starting at a large distance, integrates it outwards from a small distance to the same 
matching point, and from the difference between the logarithmic derivatives at this point 
an improved value of the energy is found. Our main innovation to this scheme is to replace 
the usual finite difference method of solving the Schrodinger equation by a method which 
solves the corresponding integral (Lippman-Schwinger) equation. That method expands 
the wave function in each radial partition in terms of Chebyshev polynomials, and solves 
matrix equations for the coefficients of the expansion. Increased accuracy is obtained by 
this procedure for three reasons: a) the solution of an integral equation is inherently more 
accurate than the solution of a differential equation; b) by using integral equations, the 
derivatives of the wave function required at the internal matching points can be expressed 
in terms of integrals that are more accurate than calculating the derivatives by a numerical 
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three- or five point formula, and c) because of the spectral nature of the expansion of the 
wave function in each partition, the length of each partition can be automatically adjusted 
in order to maintain a prescribed accuracy. This last property enables the S-IEM to treat 
accurately the abrupt transition of the wave function from the repulsive core region into 
the attractive valley region. This feature, once applied to the solution of the three-body 
problem, is also of importance in the exploration of the Efimov states [iG^ . 

To illustrate this method, the binding energy of the He dimer has been calculated, based 
on the TTY potential given by Tang, Toennies, and Yiu The result is close to the ones 
quoted in the literature, as displayed in Table IVTTI Additional numerical properties of the 
S-IEM have been explored by means of the He — He example. The accuracy of the binding 
energy was found to faithfully track the input value of the tolerance parameter, as is shown 
in Table IIIll The meshpoint economy of the method is very good. For an accuracy of three 
significant figures, the number of meshpoints needed in the radial interval between and 
3, 000 a.u. required only 208 mesh points. After an addition of 70 meshpoints, the accuracy 
increased to six significant figures. 

One of the authors (GR) acknowledges useful conversations with F. A. Gianturco, W. 
Glockle, 1. Simbotin, W. C. Stwalley, and K. T. Tang 
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Appendix 1: Recursion Relations for the coefficients 
aand b. 



The recursion relation between coefficients a and b , from one partition to a neighbouring 
partition is given by Eq. (jl3j) in the text. The corresponding matrices uJi and 7^ are given 



by 







1 



(19) 



1 - {9y)^ 



{fy)^ 1 



(20) 



1 







where 




(21) 



16 



{f^)i = / fiir)V{r)zi{r)dr 
'ti-i 

u 

{gy)i = I gi{r)V{r)yi{r)dr 
'ti-i 



gi{r)V{r)zi{r)dr. 



(22) 



(23) 



(24) 



Equation (|T!^ enables one to march outward 



i+l/ 



Ei+i/Ei 
Ei/Ei+i 







li 




bo,i 



or inward 



The integration outward is started at the innermost partition i = 1 with 



ai,i 


= 7-^ 


bi,i 





Ei/Ei+i 










Ei+i/Ei 




bi,i+i 



(25) 



(26) 



ao,i 




1/^1 


bo,i 








(27) 



and the integration inwards is started at the outermost partition (ending at T), for which 
the coefficients am and bm are given as 



(28) 



If the calculation of positive energy Sturm-Liouville functions is envisaged, whose asymptotic 
behavior is exp(2A;r) and approach for r — 0, then ai^m = i and bj^m = 1, while ao,i = 1 
and bo^i = 

The values of the functions y and z and their derivatives at upper and lower end-points 
ti and ti_i of partition i, required in the evaluation^ of Eq. ((Tj), are obtained from integral 
equations that these functions obey. The result is 










bl,m 




Em 



(29) 



z,{U)=g,{U){l-{fz),), 
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(30) 



yi{ti-i) = MU-i){l - {gy)i), 



(31) 



ZiiU^i) = giiti-i) - fiiti^i) {gz)^ . (32) 

Expressions for the derivatives of y and z at upper and lower end-points ti and of 
partition i are obtained by replacing functions / and g by their respective derivatives in 
the above equations. Since derivatives of the functions / and g are given analytically, the 
values of the derivatives of y and z at the end-points are obtained without loss of accuracy, 
contrary to what is the case when finite difference methods are employed 

Appendix 2: Units 

The transition from Eq. to the dimensionless Eq. (0) is accomplished by transforming 
the potential and the energy into dimensionless quantities according to Eqs. (fT^ and (fT3j). 
The normalization constant is given by 

^ 2M 2 ^ 2M 

Q = -T^ao X 2M = , (33 

where oq is the Bohr radius, 2R is the atomic energy unit (M ~13.606el^), h is Plank's 
constant divided by 27r, M is the reduced mass of the colliding atoms , and rrie is the mass 
of the electron. 

For the case of two colliding He atoms interacting via the TTY potential we take the 
mass of the He atom as given in Ref. i.e., h'^/Miue = 12.12 K for which the value 
of Q is 

Q = 7295.8356 (34) 

Once is obtained as the eigenvalue of equation Q, then the corresponding value of E in 
units of eV is given by 

E = -— X (27.211396) eV (35) 

It is also useful to express the energy in units of the Boltzman constant, denoted by K in 
atomic language. In this case E is given as 

- k'^ 27.211396 , , 

E = -— X K. (36) 

Q 8.617385 x IQ-^ ^ ^ 
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s 




Iters (from ((T))) 





3.0 E -3 


-2.5002592843 E-3 


1 


5.5002592823 E-3 


-1.0967998971 E-5 


2 


5.5112272813 E-3 


-2.0105203008 E - 10 


3 


5.5112274823 E-3 


-4.9700035857 E - 16 



TABLE VII: Convergence of the iterations for the wave number . The quantitie after the letter E 
denote the powers of 10 by which the quantities are to be multiphed. 



Rcutiao) 


B.E.{m.K) 


<r> (ao) 


2.0 


1.3146101 


97.7419 


2.5 


1.3146101 


97.7419 


3.0 


1.3146143 


97.7418 


3.5 


1.3219315 


97.4935 



TABLE VIII: Sensitivity of the He-He Binding Energy to the value of the cuting-off radius Rcut 

Appendix 3: Accuracy Considerations 

The quantities required for Eqs. ()25|1 and ()2()|1 are known to the same accuracy as 
the functions y and z in each partitions, given by the value of the tolerance parameter e. 
The propagation of the coefficients aj and hi across the partitions involves as many matrix 
inversions and multiplications in Eqs. ()25|) and OUj) as there are partitions, and thus the 
accuracy of Hs for each iteration, given by Eq. ((7j), is reduced by toZx number of partitions. 
The number of partitions is approximately 30, hence for e = 10^^^ the accuracy of the final 
wave number eigenvalue k is expected to be better than 10"^''. 

The rate of convergence of the iterations is shown in Table IVIII . The sensitivity of 
the binding energy to the values of Rcut is given in Table lyTTTl . The table shows that the 
repulsive core has a non- negligible effect in the 7th significant figure beyond 2.5 a^. 
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